import sys
import cplex
import numpy as np
from cplex.exceptions import CplexError

ca = cplex.Cplex()
ca.read(sys.argv[1])
#ca.read("bell5.mps")

#copy ca to c for self processing
c = cplex.Cplex(ca)

# change the type from MILP into LP (0:LP, 1:MILP)
c.set_problem_type(type=0)

# c.solve()
# d.solve()

# three vectors to track the progress of tree search process
# use three list to store the information
# P: path vector
# S: sign vector (the sign indicate the bound value)
P = []
S = []
s = 0
# declare global z and node z
z_global = -10000000
z_node = 0

# c.variables.set_lower_bounds to fix values at 0 or 1
# c.solution.get_status() to know the status for logical control  
# 1. optimal
# 2. unbounded
# 3. infeasible 
#
#

# convert output from cplex into numpy ndarray format
#       x = np.array(x)      
#
# Branching variable selection:
# For branching on the most fractional variable
# Test whether solution of LP is integer feasible
# if yes, fathom the node by integrality
# if false, we should trace the index of most fractional
# variable
#
#      s = argmax(np.minimum( x - np.floor(x) , np.ceil(x) - x) )
#
# then branch up by modifying upper bound of x_s
#
#      c.variables.set_upper_bounds(s, 1)
#
# or branch down by modifying lower bound of x_s
#
#      c.variables.set_lower_bounds(s, 0)

def backtrack(c, P, S):
    #print len(P), len(S)
    #print P[-1], S[-1]

    if len(P) == 0:
        return []

    # check the sign last elements of P
    while P[-1] > 0:
        # skip through all fathomed nodes
        if S[-1] == -1:
            c.variables.set_upper_bounds( abs(P[-1]) - 1, 1 )
        else:
            c.variables.set_lower_bounds( abs(P[-1]) - 1, 0 )
        P.pop()
        S.pop()
        if len(P) == 0:
            break

    # if back to the root, return 
    if len(P) == 0:
        return []
    
    # then change the sign of S
    if S[-1] == -1:
        c.variables.set_upper_bounds( abs(P[-1]) - 1, 1 )
        c.variables.set_lower_bounds( abs(P[-1]) - 1, 1 )
    else:
        c.variables.set_lower_bounds( abs(P[-1]) - 1, 0 )
        c.variables.set_upper_bounds( abs(P[-1]) - 1, 0 )
        
    P[-1] = -P[-1]
    S[-1] = -S[-1]

    print 'P', P, 'S', S
    
    return 
#End of backtrack function

iter = 0

while True:
    # be ready to solve the node LP
    iter = iter + 1

    c.solve()
    print 'Iteratation', iter
    print 'status', c.solution.get_status_string(), c.solution.get_status()
    # check the solution status
    if c.solution.get_status() == 3:
        # fathom by infeasibility
        # update P and S and add modify c
        backtrack(c, P, S)
        print 'visit here'
    else:
        # get solution at node LP
        x = np.array(c.solution.get_values())

        # check integer feasibility
        if (x - np.floor(x)).any():
            # integer infeasible then check the bounds
            z_node = c.solution.get_objective_value()

            if z_node <= z_global:
                # fathom by bounds and then backtrack
                backtrack(c, P, S)                
            
            # then choose the branching variable
            s = np.argmax(np.minimum( x - np.floor(x) , np.ceil(x) - x) )
            ### update P and S and modify c and go the next node always go left
            P.append(-s - 1)
            c.variables.get_lower_bounds()
            print 's', s
            # go left
            # S.append(-1)
            # c.variables.set_upper_bounds(s, 0)
            #go right
            S.append(1)
            c.variables.set_lower_bounds(s, 1)
        else:
            # integer feasible then update the node optimal value
            z_node = c.solution.get_objective_value()
            z_global = max(z_node, z_global)

            # fathom by integrality
            backtrack(c, P, S)
    print P
    print S
    print x
    print 'z_node', z_node
    print 'z_global', z_global

    if len(P) == 0:
        break

# post processing    
#sol = c.solution

#print   
# print out problem info
#print c.linear_constraints.get_rows()

# solution.get_status() returns an integer code
#print "Solution status = " , sol.get_status(), ":",
#print sol.status[sol.get_status()]
#print "Solutions: ", sol.get_values()
#print "Objective values: " , sol.get_objective_value()

